 function [GG,impact,C,eu,SDX,ZZ,initss,ssvec,flag,ssnames,stanames,shonames]=modBrookingBig(parest,~,~,calval,posest,poscal)
% =========================================================================
% This is DtFacInfExpFameCRGDP 
% 
%% *I. Description for CHIMOD DT FACP X SPREAD NEWS4* 
% 
% Model code for Chicago Fed model (CHIMOD) with 
% Double Trend (DT) Factor structure for prices (FACP) 
% Separates the ARMA(1,1) wage markup/labor disutility shock into a 
% white noise wage markup and an AR(1) labor disutility shock (X) 
% with a spread observable (SPREAD) and allowing for policy news shocks
% upto 4 periods ahead (NEWS4). 
% 
%% *II. Input* 
%
% There are either 1 or 4 inputs 
% *CASE 1* If there is 1 input, then PAREST should be the whole parameter
% vector. 
% 
% *CASE 2* If there are 4 inputs, the parameter vector PARAM is created
% from the estimated and calibrated values. 
% PAREST are the estimated values in rows POSEST and CALVAL are the
% calibrated values in position POSCAL. 
%
%% *III. Output*  
%
% The model to be solved is written, inside this code, as $$ A_{0}Y_{t+1} +
% A_{1}Y_{t}=A_{2}Y_{t-1} + PSI \eta_{t} $$  where the A matrices are
% filled with the linearized equation and entries depend on the parameter
% vector PARAM. 
%
% The solution of the model is a state space system characterized by 
%
% *Transition Equation* $$ S_{t} = GG S_{t-1} + RR \eta_{t} $$ and variance covariance matrix, 
% where $$ S_{t} $$ is the vector of all
% variables (endogenous and exogenous) of the model of dimension [NY 1],
% and $$ \eta $$  is the vector of exogenous innovations of dimension [NX 1] with 
% variance covariance matrix $$ V(\eta) = SDX'SDX $$. 
% Matrix dimensions: 
% GG [NY NY], IMPACT [NY NX], SDX [NX NX].  
%
% *Measurement Equation* $$ Y_{t} = ZZ S_{t} + ZZ C  $$ , 
% where $$ Y_{t} $$ is the data vector of dimension [NZ 1], $$ ZZ $$ is a selection matrix indicating 
% which states are observed, and $$ CC $$ is a vector of constants of
% dimension [NY 1]. 
% 
% EU is a [2 1] vector equal to [1;1] if there is a unique stable RE
% solution 
% 
% NETA GEV and COF are obsolete outputs that are still called in some codes
% and will be deleted in the new streamlined version of the estimation
% routine. 
%
%% *IV Notes on Multiple Prices* 
% 
% This version accomodates as observables an Investment Deflator and an additional NSERP prices. 
% Of these NSERP prices the first one should be the GDP deflator and the
% remaining NSERP-1 series multiple observables for the Consumption
% Deflator (e.g. chain weighted model consistent, Core PCE, Core CPI). 
% 
%% *V. Notes on Parametrization* 
% 1. Taylor rule written in terms of annualized output and inflation
% growth, with an inflation drift. 
% 2. Wage markup and labor disutility shock split. 
% 3. The composite steady state (SS) growth rate of the trend and the (SS)
% growth rate of ISTS progress are controlled by separate parameterss
% *gamma100* and *gammamiu100* rather than two parameters connected with
% ALPHA. This facilitates putting priors. 
% 4. AR coefficients of idiosyncratic errors are given by 2*PAR-1, where
% PAR is the paramater plugged into solution. Hence while PAR should be in
% the [0,1] interval, the AR coefficients are allowed to go negative. 
% 
%% *IV.a Model and data means* 
% 
% The model based means in the observation equation are given by the
% parameters *gammamiu100* and *pss100* which correspond to the mean growth
% rate of the relative price of investment to consumption and the mean
% growth rate of the consumption deflator. In addition, we allow the means
% of the observed price series to differ from these model based means
% according to extra parameters, intended to capture the discrepancy
% between model and data. 
% 
% Hence the model based mean of the Investment deflator is obtained as
% *pss100-gammamiu100-CMIU*.
% Similarly, the constants for the Consumption prices  is given by
% *pss100+CP* where CP differs by price. For simplicity this is also the
% constant for the GDP deflator, rather than aggregating that for
% investment and consumption; historically for instance Core PCE and the
% GDP deflator have grown at very similar rates. 
% 
%% *V. Related models* 
% 
% *MOD CHIMOD DT FAC P3 X SPREAD NEWS4 _ OBS:* is used to declare
% additional variables not used in estimation but annalized (e.g. output
% gap, real rate). 
% 
% *MOD CHIMOD DT FAC P3 SPREAD NEWS4:* very similar code but without split
% of wage markup and labor distulity shocks 
% 
%% *Version*
% Created by Alejandro Justiniano 12/13/2010  
% 
% Last Modified 12/13/2010. Run code ZB_CHECKWAGESPC to ensure alternative
% form of the wage Phillips curve matches model without split of wage markup and labor
% disutility shock. Also IRF of I to neutral shock is reasonable. 
% =========================================================================

if nargout < 11 
    flag_repstanames=0; 
else 
    flag_repstanames=1; 
end 
if nargout < 10 
    flag_repssnames=0; 
else
    flag_repssnames=1; 
end

flag.ssok=1;
flag.solver=1; 

initss=[]; 
ssnames=[]; 
ssvec  =[]; 


% if nargin~=1 && nargin~=4 
%     error('Need either 1 or 4 inputs') 
% end 

% =========================================================================
%% *PART 1. Declare All Variables in the Model Solution $$ S_{t} $$* 

%__________________________________________________________________________
% 1.a Endogenous Variables 
y=1;            % output
k=2;            % capital
L=3;            % hours
Rk=4;           % return on capital
w=5;            % real wages
x=6;            % marginal utility of labor
p=7;            % inflation
s=8;            % marginal cost
lambda=9;       % multiplier
c=10;           % consumption
R=11;           % interst rate
u=12;           % capital utilization
%% Delete PHI replace with q
q=13;           % q: real price of installed capital in C units
i=14;           % investment
kbar=15;        % "gross" capital
%__________________________________________________________________________
% 1.b Exogenous driving processes
mp=16;
z =17;          % Productivity shock
g =18;          % Public spending shock
miu=19;         % Relative price (non stationary), this is Growth Rate
lambdap=20;     % Markup shock AR(1) 
psi=21;         % Preference shock AR(1) 
b=22;           % Discount factor shock 
upsilon=23;     % MEI, stationary, investment shock
pitarg=24;      % Inflation Drift
wmarkwn=25;     % White Noise Wage Markup 
%__________________________________________________________________________
% 1.c Observable growth rates (nominal) 
ynomobs=26; 
cnomobs=27;
inomobs=28; 
wnomobs=29; 
%__________________________________________________________________________
% 1.c.a Auxiliary Variables or Observable Prices
yrobs   =30; 
gdp     =31; 
dpinv   =32; 
gdpdef  =33; 
%__________________________________________________________________________
% 1.d Lags
yrobs_1=34;   % Real Observable Output (adjusted for utilization) lag1      
yrobs_2=35;   % Real Observable Output (adjusted for utilization) lag2
yobsan =36;   % Annualize Observable Output 
p_1    =37; 
p_2    =38; 
dpma   =39;   % Four quarter inflation  
% New variables with the FA
nw       =40;    % Real Net Worth  
rfa      =41;    % Return on Capital, do not confuse with Rk, rental rate of capital
rfaEtpone=42;    % E[Return on Capital(t+1)|t ]
varsig   =43;    % Shock to net Worth 
b2y      =44;    % Borrowing to GDP
db2y     =45;    % First Difference Borrowing to GDP
dagtrend =46;  % First difference of the aggregate trend 
spread   =47; 

R_lead1=48;  
R_lead2=49; 
R_lead3=50; 
R_lead4=51;
%__________________________________________________________________________
% 1.e 1 through 4 step ahead expectations of FFR 
startstar=51;
% 1.f Policy signals
stsig=startstar; 
sigR_0=stsig+1; 
sigR_1=stsig+2; 
sigR_2=stsig+3; 
sigR_3=stsig+4; 
sigR_4=stsig+5; 

NY=stsig+5; 


%% Removed when shutting down expected inflation 
% % 1.g Inflation Expectations 
% posExpInf=stsig+5+1; 
% %% Expectations of inflation 1 through 40 quarters ahead
% pe_vec=(posExpInf):(posExpInf+39); 
% expectedPi40=posExpInf+40; 
% NY=posExpInf+40;

% =========================================================================

% =========================================================================
%% *PART 2. Declare Innovations to the exogenous variables $$ \eta_{t} $$* 

nsignals = 4;    % Number of policy signals 
NX=11+ nsignals; 
%__________________________________________________________________________
% 2.a Structural Disturbances 
Rs=1;           % MP
zs=2;           % Technology 
gs=3;           % Public spending
mius=4;         % Investment     
lambdaps=5;     % Mark-up
psis=6;         % Leisure preference
bs=7;           % Discount factor
upsilons=8;     % External Finance Premium 
shpitarg=9;     % Inflation Drift
wmarkwnsh=10;   % Wage markup 
varsigsh =11;   % Shock to net Worth
%__________________________________________________________________________
% 2.b Policy Signals 
%__________________________________________________________________________
% 2.b Policy Signals 
Rs_lag1=12;    % 1 period ahead
Rs_lag2=13;    % 2 periods ahead
Rs_lag3=14;    % 3 periods ahead
Rs_lag4=15;    % 4 periods ahead
% =========================================================================

% =========================================================================
%% *PART 3. Declare Parameters* 
ncof=        74;    % Number of coefficients 
numpar=      ncof;  % Equal to Parameters, distinction irrelevant here 
if nargin < 4 || isempty( calval ) 
    param=parest  ; 
else
    param=zeros(numpar,1); 
    param(posest)=parest; 
    param(poscal)=calval; 
end; 
%__________________________________________________________________________
% 3.1 Parameters for all but multiple prices, GDP and I deflators 

alpha=param(1);         % share of capital in the prod. function
delta=param(2);         % capital depreciation rate
iotap=param(3);         % price indexation (=0 is static indexation, =1 is dynamic)
iotaw=param(4);         % wages indexation
gammastar100=param(5);  % steady state growth rate of technology
gammamiu100=param(6);   % steady state growth rate of capital embodied technology
h=param(7);             % habit formation parameter
lambdapss=param(8);     % steady state of mark-up shock (pins down steady state of wages)
lambdaw=param(9);       % wage mark-up
Lss=param(10);          % steady state for leisure (the ss for leisure is pinned down by psiss. convenient to parameterize the model in terms of Lss)
pss100=param(11);       % steady state quarterly inflation (multiplied by 100)
Fbeta=param(12);        % weird parameter of SW that is a function of beta and the steady state quarterly real rate of interest (ensures that beta is less than 1)
gss=param(13);          % steady state of the public spending shock
niu=param(14);          % curvature of the utility function for leisure
xip=param(15);          % price stickiness
xiw=param(16);          % wage stickiness
chi=param(17);          % elasticity of the capital utilization cost function
S=param(18);            % investment adjustment cost
fp=param(19);           % reaction to annualized inflation
fy=param(20);           % reaction to annualized output growth 
rhoR=param(21);         % policy smoothing  
% AR coefficients of exogenous prorcesses 
rhoz=param(22); 
rhog=param(23); 
rhomiu=param(24); 
rholambdap=param(25); 
rhopsi=param(26); 
rhob=param(27); 
rhopitarg=param(28); 
rhomp=param(29); 
rhoupsilon=param(30); 
%% New Parameters with the FA
rhovarsig =param(31);   % Persistence shock to net worth
BNratio   =param(32);   % Capital to NW ratio in SS 
KNratio   =BNratio+1; 
FKN100    =param(33);   % 100*log( F(K/N) ), approximately the spread    
tau       =param(34);   % Elasticity of external premium w.r.t. KQ/N
zeta      =param(35);   % Survival Probability of entrepeneurs 
% Std of structural shocks 
sdstruct=param(36:46); %% 11 Shocks, last one is Net Worth VARSIG 
stdSpreadME= param(47); %% Measurement error spread
stdB2YME= param(48); 
% Std Policy Signals
stdsignals=1e-7*ones(nsignals,1); 
stdsignals(1:4)=param(49:52); 

% stdsignals2 is for the counterfactuals

%% Parameters for inflation expectations 
posCPI=param(53); 
expInfLoading=param(54); 
end_nonp=54; 
% _________________________________________________________________________
% 3.2  Begin Parameters for Multiple Prices 
%__________________________________________________________________________
% 3.2.1  How to aggregate gdp deflator
flag_wadd1  =param(end_nonp+1); 
nserp       =param(end_nonp+2); 
% __________________________________________________________________
% 3.2.1 Declare positions for 
% Rows for 
% a) loadings 
% b) ar ID errors 
% c) constants 
% d) volatilities 
%
% *Note:* There are NSERP+1 loadings. First 2 are for the GDP Deflator.
% 1st of W(1)*Price consumption.2nd on W(2)*Prince Investment. 
% Remaining are for other prices on the Consumption Deflators 
row_load=(end_nonp+3):(end_nonp+nserp+3); 
row_ar  =row_load(end)*ones(1,nserp)+(1:nserp); 
row_std =row_ar(end)*ones(1,nserp)+(1:nserp);
row_c   =row_std(end)*ones(1,nserp)+(1:nserp); 

%__________________________________________________________________________
% 3.2.3 vector of factor loadings 
bet_gdpdef  =param(row_load(1:2)); 
bet_c       =param(row_load(3:end)); 

bet_CPI     =bet_c(posCPI-1); 
if abs(expInfLoading-999) < 0.01; 
    expInfLoading=bet_CPI;
end



%__________________________________________________________________________
% 3.2.4 vector of idiosyncratic errors 
rhoid=2*param(row_ar)-1; 
if any(rhoid < -0.995)~=0 || any(rhoid > 0.995)~= 0 
    disp('Explosive AR ID roots'); 
    GG=[];C=[];RR=[];eu=[0;0];SDX=[];ZZ=[];
    return
end

%__________________________________________________________________________
% 3.2.5 vector of constants 
cprices =param(row_c); 
if all( cprices ~= 0 ) 
    error('At least one of the constants must be set equal to 0')
end
constantCPI=cprices(posCPI); 

%__________________________________________________________________________
% 3.2.6 vector of idiosyncratic errors
sdid             =param(row_std); 

%__________________________________________________________________________
% 3.2.7 constants for I deflator 
cons_idef        =param(row_c(end)+1:end);
if length(cons_idef)> 1 
    error('Wrong Dimension of Investment Deflator') 
end 

%==========================================================================
%% *PART 4. Compute Steady State (SS)* 
 
gammastar=gammastar100/100;
gammamiu=gammamiu100/100;
beta=100/(Fbeta+100);
rss=exp(gammastar)/beta-1;
rss100=rss*100;
gss=1/(1-gss);
%Rkss=exp(gammastar+gammamiu)/beta-1+delta;

%% Additional variables from the FA
% SS return on Capital in FA model 
rfaSS=exp(FKN100/100)*(1+rss)-1;
nomss=rss+pss100/100; 
% SS spread 
spreadSS=rfaSS-rss;
% Rental rate of capital 
Rkss=exp(gammamiu)*(1+rfaSS)-(1-delta); 

sss=1/(1+lambdapss);
wss=(sss*((1-alpha)^(1-alpha))/((alpha^(-alpha))*Rkss^alpha))^(1/(1-alpha));
kLss=(wss/Rkss)*alpha/(1-alpha);
FLss=(kLss^alpha-Rkss*kLss-wss);
yLss=kLss^alpha-FLss;

% Scale of the economy is fixed at 1, i.e. independent of LSS which now
% only enters observation equation for hours 
Lsscale=1;
kss=kLss*Lsscale;
iss=(1-(1-delta)*exp(-gammastar-gammamiu))*kss*exp(gammastar+gammamiu);
F=FLss*Lsscale;
yss=yLss*Lsscale;
css=yss/gss-iss;

% SS net worth 
nss=kss/KNratio; 
% SS entrepeneur transfer
xss=nss*(1-exp(-gammastar)*zeta*(1+rfaSS))/(1-zeta); 

ssvec  =[css/yss;iss/yss;kss/yss;rss100;100*Rkss;100*rfaSS;100*spreadSS;kss/nss]; 
if  flag_repssnames==1;
    ssnames={'c/y'  ;'i/y'  ;'k/y';  'rss100';'Rkss100';'RkFAss100';'100*Spread';'k/n'};
    %printcell([ssnames num2cprec(ssvec)]); 
end
    
%%==========================================================================

%==========================================================================
%% *PART 5. Declare System Matrices* 
% Model to be solved is written as 
% $$ A_{0}Y_{t+1} + A_{1}Y_{t}=A_{2}Y_{t-1} + PSI \eta_{t} $$  
A0  = zeros(NY,NY) ;
A1  = zeros(NY,NY) ;
A2  = zeros(NY,NY) ;
PSI = zeros(NY,NX) ;
%==========================================================================
expg=exp(gammastar);

%% *PART 6. Begin Filling Equation Rows* 

%% Eq 1, Output (y) from production function 
A1(y,y)=1;
A1(y,k)=-((yss+F)/yss)*alpha;
A1(y,L)=-((yss+F)/yss)*(1-alpha);

%% Eq. 2, Effective Capital (K) 
A1(k,k)=1;
A1(k,u)=-1;
A1(k,z)=1;
A2(k,kbar)=1;
A1(k,miu)=alpha/(1-alpha)+1;

%% Eq. 3 Hours (L), from cost minimization 
A1(L,Rk)=1;
A1(L,w)=-1;
A1(L,k)=1;
A1(L,L)=-1;

%% Eq. 4 Return to Capital (Rk) from utilization FOC
A1(Rk,Rk)=1;
A1(Rk,u)=-chi;

%% Eq. 5, Real Wages (w) PC (with split of shocks)
slopewages=(1-beta*xiw)*(1-xiw)/((1+beta)*xiw*(1+niu*(1+1/lambdaw)));
A1(w,w)=1;
A0(w,w)= -beta/(1+beta);
A0(w,p)= -beta/(1+beta);
A1(w,x)=slopewages;
A1(w,z)=(1+iotaw*beta-rhoz*beta)/(1+beta);
A1(w,miu) =(alpha/(1-alpha))*((1+iotaw*beta-rhomiu*beta)) *(1/(1+beta)); 

A1(w,p)=(1+iotaw*beta)/(1+beta);
A2(w,p)=iotaw/(1+beta);
A2(w,w)=1/(1+beta);

A2(w,z)=iotaw/(1+beta);
A2(w,miu)=(alpha*iotaw)/( (1+beta)*(1-alpha) ); 

A1(w,wmarkwn)= -1; 


%% Eq. 6, Wage Gap (x), (MRS-Real Wage)
A1(x,x)=1;
A1(x,w)=-1;
A1(x,lambda)=-1;
A1(x,b)=1/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));
A1(x,L)=niu;
A1(x,psi)=1; 

%% Eq 7, Prices (p) PC 
A1(p,p)=1;
A0(p,p)=-beta/(1+iotap*beta);
A2(p,p)=iotap/(1+iotap*beta);
A1(p,s)=-((1-beta*xip)*(1-xip)/((1+iotap*beta)*xip));
A1(p,lambdap)=-1;       

%% Eq 8, Marginal Cost (s)
A1(s,s)=1;
A1(s,Rk)=-alpha;
A1(s,w)=-(1-alpha);

%% Eq 9, Lagrange Multiplier IBC (Lambda) 
A1(lambda,lambda)=1;
A1(lambda,R)=-1;
A0(lambda,lambda)=-1;
A0(lambda,p)=1;
A1(lambda,z)=rhoz;
A1(lambda,miu)=(rhomiu)*alpha/(1-alpha);

%% Eq 10, Consumption (c)
A1(c,lambda)=(expg-h*beta)*(expg-h);
A1(c,b)=-(expg-h*beta*rhob)*(expg-h)/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)); 
A1(c,z)=-(beta*h*expg*rhoz-h*expg);
A1(c,miu)=(-(beta*h*expg*rhomiu-h*expg))*alpha/(1-alpha);
A1(c,c)=(expg^2+beta*h^2);
A0(c,c)=-beta*h*expg;
A2(c,c)=expg*h;

%% Eq 11, Taylor Rule for Nominal Interest Rate (R)
A1(R,R)=1;
% 10.1 Systematic Component 
A2(R,R)=rhoR;
A1(R,dpma)  = -(1-rhoR)*fp; % Annualized Inflation 
A1(R,yobsan)= -(1-rhoR)*fy; % Annualized Observable Output 
% 10.2 Non-Systematic Component 
A1(R,mp)      = -1;         % Contemporaneous Innovation 
A1(R,sigR_0)  = -1;         % Signals 
A1(R,pitarg)  =(1-rhoR)*fp; % Inflation Drift

%% Eq. 12, Utilization (U) from market clearing condition 
A1(u,c)=css/yss;
A1(u,i)=iss/yss;
A1(u,y)=-1/gss;
A1(u,g)=1/gss;
A1(u,u)=kss*Rkss/yss;

%% Eq. 13, Price of Installed capital in C units(Q) modified by FA
% Replace this equation with the FOC for the return on capital 
A1(q,rfaEtpone)= -1; 
A0(q,Rk)=Rkss/(Rkss + (1-delta) ); 
A0(q,q )=(1-delta)/(Rkss + (1-delta) ); 
A1(q,miu)= -rhomiu; 

A1(q,q)= -1; 

% Old equation 
% A1(phi,phi)=1;
% A0(phi,phi)=-beta*exp(-gammastar-gammamiu)*(1-delta);
% A1(phi,z)=rhoz;
% A1(phi,miu)=(rhomiu)*(alpha/(1-alpha)+1);
% A0(phi,lambda)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));
% A0(phi,Rk)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));


%% Eq. 14 Real return on capital expected next period  NEW 

A1(rfa,rfa)=  -1; 
A1(rfa,Rk)=Rkss/(Rkss + (1-delta) );  
A1(rfa,q)=(1-delta)/(Rkss + (1-delta) ); 
A1(rfa,miu)= -1; 
A2(rfa,q)=1;  

%% Eq. 15, Definition of the external finance premium NEW
% Note that this combined with the two previous equations defines 
A1(rfaEtpone,rfaEtpone)= -1; 
A0(rfaEtpone,p  )= -1; 

A1(rfaEtpone,R)=1;  
A1(rfaEtpone,kbar)=tau; 
A1(rfaEtpone,q)=tau; 
A1(rfaEtpone,nw)= -tau; 
A1(rfaEtpone,upsilon)=1; 

A1(spread,spread)= -1; 
A1(spread,rfaEtpone)=1; 
A1(spread,R)= -1;  
A0(spread,p  )=1; 

%% Eq. 16 Net Worth NEW
A1(nw,nw    )=exp(gammamiu)/( zeta*(1+rfaSS) ); 
A1(nw,rfa   )= -KNratio; 
A1(nw,varsig)= -1/A1(nw,nw); 

A2(nw,nw)=1; 
A2(nw,rfaEtpone)=(1-KNratio); 
A2(nw,z)= -1; 
A2(nw,miu)= -1; 

%% Debt to GDP Ratio NEW
A1(b2y,b2y)= -1; 
A1(b2y,q  )=1/(1-(1/KNratio)); 
A1(b2y,kbar)=1/(1-(1/KNratio)); 
A1(b2y,nw  )=-( 1/(KNratio-1)); 
A1(b2y,gdp)= -1; 

%% First Difference of Debt to GDP Ratio NEW
A1(db2y,db2y)= -1; 
A1(db2y,b2y)= 1;

A2(db2y,b2y)=1;

%% Eq. 17, Investment MODIFIED by FA 
% New Investment Equation 
expgmiu=exp(gammastar+gammamiu);
% Phi-Lambda=q 
A1(i,q)     = -1/(S*expgmiu^2); 

A1(i,miu)=((1-beta*rhomiu))*(alpha/(1-alpha)+1); %no normalization
A1(i,i)=(1+beta);
A1(i,z)=(1-beta*rhoz);
A0(i,i)=-beta;
A2(i,i)=1;


%% Eq 18, Capital accumulation (kbar)
% Note: Recall, effective capital above is KBAR*UTILIZATION 
A1(kbar,kbar)=1;
A1(kbar,i)=-(1-(1-delta)*exp(-gammastar-gammamiu));
%A1(kbar,upsilon)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A2(kbar,kbar)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,z)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,miu)=((1-delta)*exp(-gammastar-gammamiu))*(alpha/(1-alpha)+1);

%% Eq 19, MP Exogenous Policy Innovation (unexpected) 
A1(mp,mp)=1; A2(mp,mp)=rhomp; PSI(mp,Rs)=1; 

%% Eq 20, Exogenous Growth Rate in Neutral Technoloy 
A1(z,z)=1; A2(z,z)=rhoz; PSI(z,zs)=1;

%% Eq 21, Exogenous Process in G Spending 
A1(g,g)=1; A2(g,g)=rhog; PSI(g,gs)=1;

%% Eq 22, Exogenous Process in Investment Specific Tech or Relative Price C
% to I 
A1(miu,miu)=1; A2(miu,miu)=rhomiu; PSI(miu,mius)=1;

%% Eq 23, Exogenous Process in Price Markup 
A1(lambdap,lambdap)=1; A2(lambdap,lambdap)=rholambdap;PSI(lambdap,lambdaps)=1;

%% Eq 24, Exogenous Process in Labor Disutility (AR1) 
A1(psi,psi)=1; A2(psi,psi)=rhopsi; PSI(psi,psis)=1;

%% Eq 25, Exogenous Process in Intertemporal Preference
A1(b,b)=1; A2(b,b)=rhob; PSI(b,bs)=1;

%% Eq 26, Exogenous Process in External Finance Premium
A1(upsilon,upsilon)=1; A2(upsilon,upsilon)=rhoupsilon; PSI(upsilon,upsilons)=1;

%% Eq 27, Exogenous Process in Inflation Drift 
A1(pitarg,pitarg)=1; A2(pitarg,pitarg)=rhopitarg; PSI(pitarg,shpitarg)=1;

%% Eq 28, Exogenous Process in Wage Markup (white noise) 
A1(wmarkwn,wmarkwn)=1; PSI(wmarkwn,wmarkwnsh)=1; 

%% Eq 29, Exogeneous Process VARSIG in Net Worth 
A1(varsig,varsig)=1; A2(varsig,varsig)=rhovarsig; PSI(varsig,varsigsh)=1;

%% Eq. 30, Observable GDP Nominal (add P back)
A1(ynomobs,ynomobs) =  1;
A1(ynomobs,gdp)     = -1; 
A1(ynomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(ynomobs,p)       = -1; 

A2(ynomobs,gdp)     = -1; 

%% Eq. 31, Observable Consumption Nominal (add P back)
A1(cnomobs,cnomobs) =1;
A1(cnomobs,c)       = -1; 
A1(cnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(cnomobs,p)= -1; 

A2(cnomobs,c)= -1; 

%% Eq. 32 Observable Investment Nominal (add P back)
% In consumption units
A1(inomobs,inomobs)   =1;
A1(inomobs,i)         = -1; 
A1(inomobs,dagtrend)  = -1; 
% Add inflation, since nominal 
A1(inomobs,p)= -1; 

A2(inomobs,i)= -1; 

%% Eq. 33 Observable Hourly Wages  Nominal (add P back)
A1(wnomobs,wnomobs) =1;
A1(wnomobs,w)       = -1; 
A1(wnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(wnomobs,p)= -1; 

A2(wnomobs,w)= -1; 

%% Eq. 34 Observable Real GDP Growth (in C Units)
% Corrected observation equation for utilization 
% $$ yrobs_{t} = gdp_{t}-gdp_{t-1} + agtrend_{t} $$ 
% where AGTREND is the aggregate trend 
A1(yrobs,yrobs)=1; 
A1(yrobs,gdp)= -1;
A1(yrobs,dagtrend)= -1; 
A2(yrobs,gdp)= -1; 

%% Eq. 35 GDP (detrended), Real Detrended Output corrected for Utilization 
A1(gdp,gdp)= -1; 
A1(gdp,y)=1; 
A1(gdp,u)= -kss*Rkss/yss; 

%% Eq. 36 Growth Rate in Investment Deflator 
% $$ \Delta P^{I}_{t} = \Delta P^{C}_{t} - \mu_{t} $$
A1(dpinv,dpinv)= -1; 
A1(dpinv,p)    =  1; 
A1(dpinv,miu)  = -1; 

%% Eq. 37 GDP Deflator 
% $$ P^{GDP}_{t} = \beta_{1} P^{C}_{t} + \beta_{2} P^{I}_{t} $$ 
% 
% Recall there are 2 ways to define the weights 
%
% *First:* $$ {\frac{C}{C+I} , \frac{I}{C+I} } $$ 
% 
% *Second:* $$ {\frac{C}{Y}  , \frac{I}{Y}   } $$
%
% Note that in second case weights will not add up to one, which is not
% desirable, but opens up the door for adding a 3 deflator, i.e. price of G
% and NX. 
bet_gdpdef=bet_gdpdef(:); 
weights=zeros(1,2); 
if flag_wadd1==1
    weights(1)=css/(css+iss);
    weights(2)=iss/(css+iss);
else
    weights(1)=css/yss;
    weights(2)=iss/yss;
end
A1(gdpdef,gdpdef) = -1; 
A1(gdpdef,p)      = weights(1)*bet_gdpdef(1); 
A1(gdpdef,dpinv)  = weights(2)*bet_gdpdef(2); 

%% Eq. 38 Lag 1 Observable Real GDP Growth (in C Units)
A1(yrobs_1,yrobs_1)=1; A2(yrobs_1,yrobs)=1;

%% Eq. 39 Lag 2 Observable Real GDP Growth (in C Units)
A1(yrobs_2,yrobs_2)=1; A2(yrobs_2,yrobs_1)=1;

%% Eq. 40 Annualized Observable Real GDP Growth (in C Units)
% $$ \Delta y^{A}_{t}= \sum_{j=0}^{3}yrobs_{t-j} $$
A1(yobsan,yobsan)= -1;

A1(yobsan,yrobs)  =1; 
A1(yobsan,yrobs_1)=1; 
A1(yobsan,yrobs_2)=1;
A2(yobsan,yrobs_2)= -1;

%% Eq. 41 Lag 1 Inflation (C)
A1(p_1,p_1)=1; A2(p_1,p)=1;

%% Eq. 42 Lag 2 Inflation (C)
A1(p_2,p_2)=1; A2(p_2,p_1)=1;

%% Eq. 43 Annualized Inflation (C) 
% $$ \Delta P^{A,C}_{t}= \frac{ \sum_{j=0}^{3} \Delta P^{C}_{t-j} }{4} $$
A1(dpma,dpma)= -1;
A1(dpma,p)   = 1/4;
A1(dpma,p_1) = 1/4;
A1(dpma,p_2) = 1/4;
A2(dpma,p_2) = -1/4;

%% Eq. 44, Growth Rate in Aggregate Trend 
% $$ \Gamma_{t} = Z_{t} + \frac{\alpha}{1-\alpha} \mu_{t} $$ 
A1(dagtrend,dagtrend)= -1; 
A1(dagtrend,z)=1; 
A1(dagtrend,miu)=alpha/(1-alpha); 

%% Eq. 45-49 Monetary Policy Signals

%% Eq. 45. $$ \zeta^{4}_{t}       = \epsilon^{4}_{t} $$ 
A1(sigR_4,sigR_4)=1;PSI(sigR_4,Rs_lag4)=1; 

%% Eq. 46.  $$ \zeta^{3}_{t} = \zeta^{4}_{t-1} + \epsilon^{3}_{t} $$ 
A1(sigR_3,sigR_3)=1; A2(sigR_3,sigR_4)=1; PSI(sigR_3,Rs_lag3)=1; 

%% Eq. 47.  $$ \zeta^{2}_{t} = \zeta^{3}_{t-1} + \epsilon^{2}_{t} $$ 
A1(sigR_2,sigR_2)=1; A2(sigR_2,sigR_3)=1; PSI(sigR_2,Rs_lag2)=1; 

%% Eq. 48.  $$ \zeta^{1}_{t} = \zeta^{2}_{t-1} + \epsilon^{1}_{t} $$ 
A1(sigR_1,sigR_1)=1; A2(sigR_1,sigR_2)=1; PSI(sigR_1,Rs_lag1)=1; 

%% Eq. 49. $$ \zeta^{0}_{t} = \zeta^{1}_{t-1} $$, which implies that the history of
% policy signals regarding the current quarter is given by 
%
% $$ \zeta^{0}_{t} = \epsilon^{1}_{t-1} +  \epsilon^{2}_{t-2} + \epsilon^{3}_{t-3} +
% \epsilon^{4}_{t-4}   $$ 
% *Note:* The contemporaneous policy innovation appears 
% directly in the definition of the Taylor rule, in case it is AR(1). An
% alternative would be to make the current composite signal +
% contemporaneous innovation persistent. 
A1(sigR_0,sigR_0)=1;A2(sigR_0,sigR_1)=1;  

%% Added Expected FFR 1 through 4 quarters ahead 
A1(R_lead1,R_lead1)=1; A0(R_lead1,R)= -1; 
A1(R_lead2,R_lead2)=1; A0(R_lead2,R_lead1)= -1;
A1(R_lead3,R_lead3)=1; A0(R_lead3,R_lead2)= -1;
A1(R_lead4,R_lead4)=1; A0(R_lead4,R_lead3)= -1;


%% Removed when shutting down expected inflation 
% pe_vec=pe_vec(:); 
% A1(pe_vec,pe_vec)=eye(length(pe_vec)); 
% A0(pe_vec,[p;pe_vec(1:end-1)])= -eye(length(pe_vec)); 
% 
% A1(expectedPi40,expectedPi40)= -1; 
% A1(expectedPi40,pe_vec      )= (expInfLoading/length(pe_vec))*ones(1,length(pe_vec));  


%% *PART 7: Solve Model using AMASOLVE*
[GG,~,impact,~,~,~,~,eu]=amasolve(A0,A1,A2,zeros(NY,1),PSI) ;
% EU =[1;1] if unique and stable RE solution 
if ~isequal(eu,[1;1]); 
    flag.euok=0; 
    ZZ=[]; C=[]; SDX=[];
    return 
end 

flag.euok=1; 
%% *PART 8: Addition of Measurement Equations for Prices*

%% 8.1 Extend model solution by 2*NSERP rows 
% Number of rows to add to the solution 
nadd=nserp*2; 
% Part 1: Position of prices 
posp =(NY+1):(NY+nserp); 
% Part 2: Position of ID errors 
posid=(posp(end)+1):(NY+2*nserp); 
GG=blkdiag(GG,zeros(nadd)); 
% Matrix of impact coefficients 
RR=zeros(NY+nadd,NX+nserp); 
RR(1:NY,1:NX)=impact; 

%% 8.2 Idiosyncratic errors 
% Which evolve as $$ \nu^{i}_{t} = \rho  \nu^{i}_{t-1} + \varsigma^{i}_{t} $$
% as independent AR(1) processes 
GG(posid,posid   )=diag(rhoid); 
RR(posid,NX+1:end)=eye(nserp); 

%% 8.3 General Structure of Observable Price Equation 
% 
% For price $$ P^{i,C}_{t} = \beta^{i} P^{C}_{t} + \nu^{i}_{t} $$ written in terms
% of lagged variables as 
% $$ P^{i,C}_{t} = \beta^{i} G_{P^{C}} S_{t-1} + \beta^{i} R_{P^{C}}
% \eta_{t} + \rho \nu^{i}_{t-1} + \varsigma^{i}_{t} $$ ,
% where $$ \beta^{i} $$ is the series specific factor loading, and, $$ G_{P^{C}}  R_{P^{C}} $$ correspond to the rows of the solution
% matrix for the model consumption price 

%% 8.3.a First row is GDP 
GG(posp(1),1:NY)=GG(gdpdef,1:NY);
RR(posp(1),1:NX)=RR(gdpdef,1:NX); 

%% 8.3.b Remaining rows are multiple indicators for consumption prices 
% Recall Structural part of observation equation $$ \beta^{i} G_{P^{C}} S_{t-1} +
% \beta^{i} R_{P^{C}} \eta_{t}  $$ 
GG(posp(2:end),:)       =repmat( GG(p,:)    ,[nserp-1 1] ); 
RR(posp(2:end),1:NX)    =repmat( RR(p,1:NX) ,[nserp-1 1] );
% Multiply by Factor Loadings  
bet_c=bet_c(:); 
GG(posp(2:end),1:NY)=(repmat(bet_c,[1 NY])).*GG(posp(2:end),1:NY); 
RR(posp(2:end),1:NX)=(repmat(bet_c,[1 NX])).*RR(posp(2:end),1:NX);
%% 8.3.c Add idiosyncratic errors 
GG(posp,posid)   =diag(rhoid); 
RR(posp,NX+1:end)=eye(nserp); 

% Update model dimensions after adding prices 
impact=RR; 
NY=size(GG,1); 
NX=size(impact,2); 
RR=[]; 
%% *Part 9: Addition of Measurement Equations for Spread*
% 
%% 9.1 Extend model solution by 1 row 
GG    =blkdiag(GG,zeros(2));
impact=[impact zeros(NY,2);zeros(2,NX+2)];

%% 9.2 Fill in last row for the spread 
% 
% $$ sp_{t} = \kappa \upsilon_{t} + \varsigma^{sp}_{t} $$ 
posspread=NY+1;
GG(posspread,:)       =GG(spread,:); 
impact(posspread,:)   =impact(spread,:); 
impact(posspread,NX+1)=1; 


posB2Yobs=NY+2; 
GG(posB2Yobs,:)       =GG(b2y,:); 
impact(posB2Yobs,:)   =impact(b2y,:); 
impact(posB2Yobs,NX+2)=1; 

% Update model dimensions after adding spread
NY=NY+2;
NX=NX+2;

%% *PART 10: Declare Cholesky of Variance Covariance Matrix* 
% Recall $$ V = SDX'SDX $$ and the order of the shocks is 
% 1. Structural 2. Policy Signals 3. Idiosyncratic Prices 4. Measurement
% Error Spread 
SDX=diag([sdstruct(:);stdsignals(:);sdid(:);stdSpreadME;stdB2YME]);
if size(SDX,1)~=NX;error('Dimensions of errors and SDX do not match');end

%% *PART 11: Declare ZZ, Matrix of Pointers to Observables*
% Ordering 
% 1) Nominal GDP Growth. 
% 2) Nominal C   Growth.
% 3) Nominal I   Growth.
% 4) Hours. 
% 5) Nominal W   Growth. 
% 6) Feds Fund Rate. 
% 7) Investment Deflator. 
% 8) GDP Deflator. 
% 9 though NSERP-1) Indicators of C Prices. 
% Last is Spread
% dimZ=7+nserp+3; 
% Now Reduced dimension by 1 as do not observe expected inflation, added 3
% due to expected FFR 
%
ZZposAfterP=8+nserp; 
dimZ=ZZposAfterP+5; 
ZZ=zeros(dimZ,NY); 
ZZ(1,ynomobs)=1; 
ZZ(2,cnomobs)=1; 
ZZ(3,inomobs)=1; 
ZZ(4,L)=1; 
ZZ(5,wnomobs)=1; 
ZZ(6,R)=1;
ZZ(7,dpinv)=1;
ZZ(8:8+nserp-1,posp)=eye(nserp);
ZZ(ZZposAfterP,posspread)=1; 
%ZZ(dimZ-1,expectedPi40)=1; 
ZZ(ZZposAfterP+1,posB2Yobs  )=1; 
ZZ(ZZposAfterP+2,R_lead1)=1; 
ZZ(ZZposAfterP+3,R_lead2)=1; 
ZZ(ZZposAfterP+4,R_lead3)=1; 
ZZ(ZZposAfterP+5,R_lead4)=1; 


%% *PART 12: Declare C, Matrix of Constants*
C=zeros(NY,1); 
C(ynomobs)=gammastar100+pss100;
C(cnomobs)=gammastar100+pss100;
C(inomobs)=gammastar100+pss100;
C(L)=Lss;
C(wnomobs)=gammastar100+pss100;
C(R)    =pss100+rss100;
C(dpinv)=(pss100-gammamiu100)+cons_idef; 
C(posp )= pss100*ones(nserp,1)+cprices(:); 
%C(expectedPi40)=pss100+constantCPI;
C(posspread   )=100*spreadSS; 
%C(posB2Yobs  )=100*log( BNratio ); 

if flag_repstanames==0; 
    return 
end 
%________________________________________________________
%% 16.1 STANAMES: Cell with Names of ALL states 
stanames1={'yhat','khat','hours','Rental rate K','what','x','dp','Marginal Cost','lambdahat','chat',...
    'nomr','Utilization ','q','ihat','kbar','mp','z','g','miu','lambdap',...
    'psi','b','upsilon ','pitarg','wmarkwn','GDPm nominal','Cm nominal','Im nominal','Wm nominal',...  % Position 29
    'GDPm real','GDPm level','I deflator','GDPm deflator','GDPm real lag1','GDP real lag2','GDP Annualized','Inflation lag1','Inflation lag2','Inflation annualized',...
    'Net worth','Return on K','Expected return on K','Shock to net worth','Borrowing to GDP model','Chg Borrowing to GDP model','Aggregate trend',...
    'Spread model','Signal0','Signal1','Signal2',...
    'Signal3','Signal4',...
    'ExpInf1','ExpInf2','ExpInf3','ExpInf4','ExpInf5','ExpInf6','ExpInf7','ExpInf8','ExpInf9','ExpInf10',...
    'ExpInf11','ExpInf12','ExpInf13','ExpInf14','ExpInf15','ExpInf16','ExpInf17','ExpInf18','ExpInf19','ExpInf20',...
    'ExpInf21','ExpInf22','ExpInf23','ExpInf24','ExpInf25','ExpInf26','ExpInf27','ExpInf28','ExpInf29','ExpInf30',...
    'ExpInf31','ExpInf32','ExpInf33','ExpInf34','ExpInf35','ExpInf36','ExpInf37','ExpInf38','ExpInf39','ExpInf40','AvgExpInf10y',...
    'GDP deflator'}; 
if nserp==3 
    stanames2={'C deflator','PCE Core','GDP def Id','C def Id','PCE Core Id'}; 
elseif nserp==4 
    stanames2={'C deflator','PCE Core','CPI Core','GDP def Id','C def Id','PCE Core Id','CPI Core Id'}; 
else 
    disp('Nprices > 4 please add names to STANAMES') 
end 
stanames3={'Spread','Borrowing to GDP'};  %,'Spread Me','RealY','RealC','RealI','RealW','SumR'}; 
stanames=[stanames1(:);stanames2(:);stanames3(:)]; 

%% 16.2 STSHOCKS: Cell with Names of Shocks in STANAMES 
% *Note:* Technically not true for the policy signals, for which looking
% only at the cummulative signals. Hence isolate policy innovations rather
% than policy states.
stshocks={'mp','z','g','miu','lambdap','psi','b','shock spread','pitarg','wmarkwn','shock nw',...
    'Signal1','Signal2','Signal3','Signal4',...
    'GDP def Id','C def Id','PCE Core Id'};
if nserp==4
    stshocks=[stshocks,'CPI Id'];
end
stshocks=[stshocks,'Spread Me','B2Y Me'];
stshocks=stshocks(:);

%% 16.3 SHONAMES: Cell with Names of Shock Innovations 
shonames={'MP','Neutral','G+NX','ISTS','Price Markup','Labor Disutility',...
    'Discount','MEI','Inflation Drift','Wage Markup','Net Worth Shock',...
    'MP Signal 1','MP Signal 2','MP Signal 3','MP Signal 4'}; 
shonames=[shonames(:);stshocks(16:end)]; 

%% 16.4 Quick check dimensions 
if length(stanames)~=size(GG,1)
    error('Mistmatch STANAMES and NY');
end
if length(stshocks)~=size(SDX,1)
    error('Mistmatch STSHOCKS and NX');
end
if length(shonames)~=size(SDX,1)
    error('Mistmatch SHONAMES and NX');
end
